﻿// -------------------------------------------------------------------------------------------------------------------------------
// <copyright file="CrcsSparseSymmetricMatrixExtensions.cs" company="Exul">
//     Copyright (c) Exul. All rights reserved.
// </copyright>
// <author>
//     Exul
// </author>
// -------------------------------------------------------------------------------------------------------------------------------

namespace ExulLibrary.Mathematics.LinearAlgebra.Sequential.DecimalExtensions
{
    #region Usings

    using System;
    using System.Globalization;

    using ExulLibrary.Exceptions;
    using ExulLibrary.Mathematics.LinearAlgebra.Exceptions;

    #endregion

    /// <summary>
    ///     Provides extension methods for <see cref="CrcsSparseSymmetricMatrix{T}"/> with <see cref="decimal"/> type parameter.
    /// </summary>
    public static class CrcsSparseSymmetricMatrixExtensions
    {
        #region Methods

        #region public

        /// <summary>
        ///     Adds the second matrix to the specified elements of the first matrix.
        /// </summary>
        /// <param name="augend">
        ///     The augend matrix.
        /// </param>
        /// <param name="augendRowsIndexes">
        ///     The collection of zero-based indexes of rows of augend elements.
        /// </param>
        /// <param name="augendColumnsIndexes">
        ///     The collection of zero-based indexes of columns of augend elements.
        /// </param>
        /// <param name="addend">
        ///     The addend matrix.
        /// </param>
        /// <exception cref="ArgumentNullException">
        ///     <paramref name="augend"/> is <see langword="null"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="addend"/> is <see langword="null"/>.
        ///     </para>
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="augendRowsIndexes"/> is <see langword="null"/>.
        ///     </para>
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="augendColumnsIndexes"/> is <see langword="null"/>.
        ///     </para>
        /// </exception>
        /// <exception cref="ElementIndexIsOutOfRangeException">
        ///     An element of <paramref name="augendRowsIndexes"/> is negative or greater than or equal to
        ///     rows count of <paramref name="augend"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         An element of <paramref name="augendColumnsIndexes"/> is negative or greater than or equal to
        ///         columns count of <paramref name="augend"/>.
        ///     </para>
        /// </exception>
        /// <exception cref="ReadOnlyMatrixElementException">
        ///     An element to assign specified by <paramref name="augendRowsIndexes"/> and <paramref name="augendColumnsIndexes"/>
        ///     is not in the portrait of <paramref name="augend"/>.
        /// </exception>
        /// <remarks>
        ///     The number of rows of <paramref name="addend"/> to add is determined as
        ///     min(<paramref name="augendRowsIndexes"/>.Length, <paramref name="addend"/>.RowsCount).
        ///     The number of columns of <paramref name="addend"/> to add is determined as
        ///     min(<paramref name="augendColumnsIndexes"/>.Length, <paramref name="addend"/>.ColumnsCount).
        /// </remarks>
        public static void Add(
            this CrcsSparseSymmetricMatrix<decimal> augend,
            int[] augendRowsIndexes,
            int[] augendColumnsIndexes,
            DenseSymmetricMatrix<decimal> addend)
        {
#if SAFE

            if (augend == null)
            {
                throw new ArgumentNullException("augend");
            }

            if (augendRowsIndexes == null)
            {
                throw new ArgumentNullException("augendRowsIndexes");
            }

            if (augendColumnsIndexes == null)
            {
                throw new ArgumentNullException("augendColumnsIndexes");
            }

            if (addend == null)
            {
                throw new ArgumentNullException("addend");
            }

#endif

            var augendSize = augend.RowsCount;
            var begins = augend.Begins;
            var indexes = augend.Indexes;
            var augendDiagonalValues = augend.DiagonalValues;
            var augendValues = augend.TriangleValues;
            var addendValues = addend.Values;

            var addendElementIndex = 0;

            var sizeToAdd = Math.Min(addend.RowsCount, Math.Min(augendRowsIndexes.Length, augendColumnsIndexes.Length));
            for (var rowIndex = 0; rowIndex < sizeToAdd; rowIndex++)
            {
                var rowElementIndex = augendRowsIndexes[rowIndex];

#if SAFE

                if ((rowElementIndex < 0) || (augendSize <= rowElementIndex))
                {
                    var argumentName = string.Format(CultureInfo.CurrentCulture, "augendRowsIndexes[{0}]", rowIndex);
                    throw new ElementIndexIsOutOfRangeException(argumentName, rowElementIndex, 0, augendSize - 1);
                }

#endif

                var rowBegin = begins[rowElementIndex];
                var rowEnd = begins[rowElementIndex + 1];

                for (var columnIndex = 0; columnIndex <= rowIndex; columnIndex++)
                {
                    var columnElementIndex = augendColumnsIndexes[columnIndex];

#if SAFE

                    if ((rowIndex == 0) && ((columnElementIndex < 0) || (augendSize <= columnElementIndex)))
                    {
                        var argumentName = string.Format(CultureInfo.CurrentCulture, "augendColumnsIndexes[{0}]", columnIndex);
                        throw new ElementIndexIsOutOfRangeException(argumentName, columnElementIndex, 0, augendSize - 1);
                    }

                    if (rowElementIndex < columnElementIndex)
                    {
                        throw new ReadOnlyMatrixElementException("augend", rowElementIndex, columnElementIndex);
                    }

#endif

                    if (rowElementIndex == columnElementIndex)
                    {
                        augendDiagonalValues[rowElementIndex] += addendValues[addendElementIndex];
                        addendElementIndex++;
                        continue;
                    }

                    var augendElementIndex = rowBegin;
                    while ((augendElementIndex < rowEnd) && (indexes[augendElementIndex] < columnElementIndex))
                    {
                        augendElementIndex++;
                    }

                    if ((augendElementIndex < rowEnd) && (indexes[augendElementIndex] == columnElementIndex))
                    {
                        augendValues[augendElementIndex] += addendValues[addendElementIndex];
                    }
                    else
                    {
#if SAFE

                        throw new ReadOnlyMatrixElementException("augend", rowElementIndex, columnElementIndex);

#endif
                    }

                    addendElementIndex++;
                }
            }
        }

        /// <summary>
        ///     Adds the second matrix multiplied by the factor to the specified elements of the first matrix.
        /// </summary>
        /// <param name="augend">
        ///     The augend matrix.
        /// </param>
        /// <param name="augendRowsIndexes">
        ///     The collection of zero-based indexes of rows of augend elements.
        /// </param>
        /// <param name="augendColumnsIndexes">
        ///     The collection of zero-based indexes of columns of augend elements.
        /// </param>
        /// <param name="addend">
        ///     The addend matrix.
        /// </param>
        /// <param name="addendFactor">
        ///     The factor of the addend matrix.
        /// </param>
        /// <exception cref="ArgumentNullException">
        ///     <paramref name="augend"/> is <see langword="null"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="addend"/> is <see langword="null"/>.
        ///     </para>
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="augendRowsIndexes"/> is <see langword="null"/>.
        ///     </para>
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="augendColumnsIndexes"/> is <see langword="null"/>.
        ///     </para>
        /// </exception>
        /// <exception cref="ElementIndexIsOutOfRangeException">
        ///     An element of <paramref name="augendRowsIndexes"/> is negative or greater than or equal to
        ///     rows count of <paramref name="augend"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         An element of <paramref name="augendColumnsIndexes"/> is negative or greater than or equal to
        ///         columns count of <paramref name="augend"/>.
        ///     </para>
        /// </exception>
        /// <exception cref="ReadOnlyMatrixElementException">
        ///     An element to assign specified by <paramref name="augendRowsIndexes"/> and <paramref name="augendColumnsIndexes"/>
        ///     is not in the portrait of <paramref name="augend"/>.
        /// </exception>
        /// <remarks>
        ///     The number of rows of <paramref name="addend"/> to add is determined as
        ///     min(<paramref name="augendRowsIndexes"/>.Length, <paramref name="addend"/>.RowsCount).
        ///     The number of columns of <paramref name="addend"/> to add is determined as
        ///     min(<paramref name="augendColumnsIndexes"/>.Length, <paramref name="addend"/>.ColumnsCount).
        /// </remarks>
        public static void Add(
            this CrcsSparseSymmetricMatrix<decimal> augend,
            int[] augendRowsIndexes,
            int[] augendColumnsIndexes,
            DenseSymmetricMatrix<decimal> addend,
            decimal addendFactor)
        {
#if SAFE

            if (augend == null)
            {
                throw new ArgumentNullException("augend");
            }

            if (augendRowsIndexes == null)
            {
                throw new ArgumentNullException("augendRowsIndexes");
            }

            if (augendColumnsIndexes == null)
            {
                throw new ArgumentNullException("augendColumnsIndexes");
            }

            if (addend == null)
            {
                throw new ArgumentNullException("addend");
            }

#endif

            var augendSize = augend.RowsCount;
            var begins = augend.Begins;
            var indexes = augend.Indexes;
            var augendDiagonalValues = augend.DiagonalValues;
            var augendValues = augend.TriangleValues;
            var addendValues = addend.Values;

            var addendElementIndex = 0;

            var sizeToAdd = Math.Min(addend.RowsCount, Math.Min(augendRowsIndexes.Length, augendColumnsIndexes.Length));
            for (var rowIndex = 0; rowIndex < sizeToAdd; rowIndex++)
            {
                var rowElementIndex = augendRowsIndexes[rowIndex];

#if SAFE

                if ((rowElementIndex < 0) || (augendSize <= rowElementIndex))
                {
                    var argumentName = string.Format(CultureInfo.CurrentCulture, "augendRowsIndexes[{0}]", rowIndex);
                    throw new ElementIndexIsOutOfRangeException(argumentName, rowElementIndex, 0, augendSize - 1);
                }

#endif

                var rowBegin = begins[rowElementIndex];
                var rowEnd = begins[rowElementIndex + 1];

                for (var columnIndex = 0; columnIndex <= rowIndex; columnIndex++)
                {
                    var columnElementIndex = augendColumnsIndexes[columnIndex];

#if SAFE

                    if ((rowIndex == 0) && ((columnElementIndex < 0) || (augendSize <= columnElementIndex)))
                    {
                        var argumentName = string.Format(CultureInfo.CurrentCulture, "augendColumnsIndexes[{0}]", columnIndex);
                        throw new ElementIndexIsOutOfRangeException(argumentName, columnElementIndex, 0, augendSize - 1);
                    }

                    if (rowElementIndex < columnElementIndex)
                    {
                        throw new ReadOnlyMatrixElementException("augend", rowElementIndex, columnElementIndex);
                    }

#endif

                    if (rowElementIndex == columnElementIndex)
                    {
                        augendDiagonalValues[rowElementIndex] += addendFactor * addendValues[addendElementIndex];
                        addendElementIndex++;
                        continue;
                    }

                    var augendElementIndex = rowBegin;
                    while ((augendElementIndex < rowEnd) && (indexes[augendElementIndex] < columnElementIndex))
                    {
                        augendElementIndex++;
                    }

                    if ((augendElementIndex < rowEnd) && (indexes[augendElementIndex] == columnElementIndex))
                    {
                        augendValues[augendElementIndex] += addendFactor * addendValues[addendElementIndex];
                    }
                    else
                    {
#if SAFE

                        throw new ReadOnlyMatrixElementException("augend", rowElementIndex, columnElementIndex);

#endif
                    }

                    addendElementIndex++;
                }
            }
        }

        /// <summary>
        ///     Decomposes the matrix by incomplete LDLT-decomposition.
        /// </summary>
        /// <param name="matrix">
        ///     The matrix to decompose.
        /// </param>
        /// <returns>
        ///     The decomposition matrix.
        /// </returns>
        /// <exception cref="ArgumentNullException">
        ///     <paramref name="matrix"/> is <see langword="null"/>.
        /// </exception>
        /// <remarks>
        ///     The method doesn't copy the input matrix.
        /// </remarks>
        public static DecompositionMatrix<decimal> DecomposeByIncompleteLdlt(this CrcsSparseSymmetricMatrix<decimal> matrix)
        {
#if SAFE

            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

#endif

            var rowsCount = matrix.RowsCount;
            var begins = matrix.Begins;
            var indexes = matrix.Indexes;
            var triangleValues = matrix.TriangleValues;
            var diagonalValues = matrix.DiagonalValues;

            var end = begins[0];

            for (var rowColumnIndex = 1; rowColumnIndex < rowsCount; rowColumnIndex++)
            {
                var start = end;

                var sum = 0.0m;

                // calculate elements in low triangle
                var elementIndex = start;
                end = begins[rowColumnIndex + 1];
                while (elementIndex < end)
                {
                    var rowSum = 0.0m;

                    var index = indexes[elementIndex];

                    var multiplicandElementIndex = start;
                    var multiplierElementIndex = begins[index];
                    var multiplierEnd = begins[index + 1];
                    while ((multiplicandElementIndex < elementIndex) && (multiplierElementIndex < multiplierEnd))
                    {
                        if (indexes[multiplicandElementIndex] == indexes[multiplierElementIndex])
                        {
                            rowSum += triangleValues[multiplicandElementIndex] *
                                      diagonalValues[indexes[multiplicandElementIndex]] *
                                      triangleValues[multiplierElementIndex];

                            multiplicandElementIndex++;
                            multiplierElementIndex++;
                        }
                        else if (indexes[multiplicandElementIndex] < indexes[multiplierElementIndex])
                        {
                            multiplicandElementIndex++;
                        }
                        else
                        {
                            multiplierElementIndex++;
                        }
                    }

                    var elementValue = (triangleValues[elementIndex] - rowSum) / diagonalValues[index];
                    triangleValues[elementIndex] = elementValue;

                    sum += elementValue * diagonalValues[index] * elementValue;

                    elementIndex++;
                }

                diagonalValues[rowColumnIndex] -= sum;
            }

            return new DecompositionMatrix<decimal>(matrix, DecompositionAlgorithm.Ldlt);
        }

        /// <summary>
        ///     Multiplies the matrix by the dense vector and stores the result in the another vector.
        /// </summary>
        /// <param name="multiplicand">
        ///     The multiplicand matrix.
        /// </param>
        /// <param name="multiplier">
        ///     The multiplier vector.
        /// </param>
        /// <param name="result">
        ///     The vector to store the result.
        /// </param>
        /// <exception cref="ArgumentNullException">
        ///     <paramref name="multiplicand"/> is <see langword="null"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="multiplier"/> is <see langword="null"/>.
        ///     </para>
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="result"/> is <see langword="null"/>.
        ///     </para>
        /// </exception>
        /// <remarks>
        ///     The number of rows of <paramref name="multiplicand"/> to multiply is determined as
        ///     min(<paramref name="multiplicand"/>.RowsCount, <paramref name="result"/>.Length).
        ///     The number of columns of <paramref name="multiplicand"/> to multiply is determined as
        ///     min(<paramref name="multiplicand"/>.ColumnsCount, <paramref name="multiplier"/>.Length).
        /// </remarks>
        public static void Multiply(
            this CrcsSparseSymmetricMatrix<decimal> multiplicand, DenseVector<decimal> multiplier, DenseVector<decimal> result)
        {
#if SAFE

            if (multiplicand == null)
            {
                throw new ArgumentNullException("multiplicand");
            }

            if (multiplier == null)
            {
                throw new ArgumentNullException("multiplier");
            }

            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

#endif

            var begins = multiplicand.Begins;
            var indexes = multiplicand.Indexes;
            var multiplicandValues = multiplicand.TriangleValues;
            var diagonalValues = multiplicand.DiagonalValues;
            var multiplierValues = multiplier.Values;
            var resultValues = result.Values;

            var matrixSizeToMultiply = Math.Min(multiplicand.RowsCount, Math.Min(result.Length, multiplier.Length));
            for (var rowIndex = 0; rowIndex < matrixSizeToMultiply; rowIndex++)
            {
                var rowSum = 0.0m;

                var multiplierElement = multiplierValues[rowIndex];

                var elementIndex = begins[rowIndex];
                var rowEnd = begins[rowIndex + 1];
                while (elementIndex < rowEnd)
                {
                    var columnIndex = indexes[elementIndex];

                    rowSum += multiplicandValues[elementIndex] * multiplierValues[columnIndex];
                    resultValues[columnIndex] += multiplierElement * multiplicandValues[elementIndex];

                    elementIndex++;
                }

                rowSum += multiplierElement * diagonalValues[rowIndex];

                resultValues[rowIndex] = rowSum;
            }
        }

        /// <summary>
        ///     Solves a system of linear algebraic equations with the specified matrix and the right-hand part and
        ///     stores the result in the right-hand part vector.
        /// </summary>
        /// <param name="matrix">
        ///     The matrix of system.
        /// </param>
        /// <param name="rightHandPart">
        ///     The right-hand part vector as input and the result vector as output.
        /// </param>
        /// <exception cref="ArgumentNullException">
        ///     <paramref name="matrix"/> is <see lagnword="null"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="rightHandPart"/> is <see lagnword="null"/>.
        ///     </para>
        /// </exception>
        /// <remarks>
        ///     The size of system to solve is determined as minimal value from size of <paramref name="matrix"/> and
        ///     length of <paramref name="rightHandPart"/>.
        /// </remarks>
        public static void SolveWithDiagonal(this CrcsSparseSymmetricMatrix<decimal> matrix, DenseVector<decimal> rightHandPart)
        {
#if SAFE

            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (rightHandPart == null)
            {
                throw new ArgumentNullException("rightHandPart");
            }

#endif

            var matrixValues = matrix.DiagonalValues;
            var resultValues = rightHandPart.Values;

            var rowsCount = Math.Min(matrix.RowsCount, rightHandPart.Length);
            for (var rowIndex = 0; rowIndex < rowsCount; rowIndex++)
            {
                resultValues[rowIndex] /= matrixValues[rowIndex];
            }
        }

        /// <summary>
        ///     Solves a system of linear algebraic equations with the specified matrix and the right-hand part and
        ///     stores the result in the right-hand part vector.
        /// </summary>
        /// <param name="matrix">
        ///     The matrix of system.
        /// </param>
        /// <param name="rightHandPart">
        ///     The right-hand part vector as input and the result vector as output.
        /// </param>
        /// <exception cref="ArgumentNullException">
        ///     <paramref name="matrix"/> is <see lagnword="null"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="rightHandPart"/> is <see lagnword="null"/>.
        ///     </para>
        /// </exception>
        /// <remarks>
        ///     The size of system to solve is determined as minimal value from columns count of <paramref name="matrix"/> and
        ///     length of <paramref name="rightHandPart"/>.
        /// </remarks>
        public static void SolveWithHighTriangle(
            this CrcsSparseSymmetricMatrix<decimal> matrix, DenseVector<decimal> rightHandPart)
        {
#if SAFE

            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (rightHandPart == null)
            {
                throw new ArgumentNullException("rightHandPart");
            }

#endif

            var begins = matrix.Begins;
            var indexes = matrix.Indexes;
            var triangleValues = matrix.TriangleValues;
            var rightHandValues = rightHandPart.Values;

            var columnsCount = Math.Min(matrix.RowsCount, rightHandPart.Length);

            var elementIndex = begins[columnsCount] - 1;

            for (var columnIndex = columnsCount - 1; 0 < columnIndex; columnIndex--)
            {
                var rowStart = begins[columnIndex];

                var currentValue = rightHandValues[columnIndex];

                while (rowStart <= elementIndex)
                {
                    var rowIndex = indexes[elementIndex];
                    rightHandValues[rowIndex] -= triangleValues[elementIndex] * currentValue;

                    elementIndex--;
                }
            }
        }

        /// <summary>
        ///     Solves a system of linear algebraic equations with the specified matrix and the right-hand part and
        ///     stores the result in the right-hand part vector.
        /// </summary>
        /// <param name="matrix">
        ///     The matrix of system.
        /// </param>
        /// <param name="rightHandPart">
        ///     The right-hand part vector as input and the result vector as output.
        /// </param>
        /// <exception cref="ArgumentNullException">
        ///     <paramref name="matrix"/> is <see lagnword="null"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="rightHandPart"/> is <see lagnword="null"/>.
        ///     </para>
        /// </exception>
        /// <remarks>
        ///     The size of system to solve is determined as minimal value from columns count of <paramref name="matrix"/> and
        ///     length of <paramref name="rightHandPart"/>.
        /// </remarks>
        public static void SolveWithHighTriangleAndDiagonal(
            this CrcsSparseSymmetricMatrix<decimal> matrix, DenseVector<decimal> rightHandPart)
        {
#if SAFE

            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (rightHandPart == null)
            {
                throw new ArgumentNullException("rightHandPart");
            }

#endif

            var begins = matrix.Begins;
            var indexes = matrix.Indexes;
            var triangleValues = matrix.TriangleValues;
            var diagonalValues = matrix.DiagonalValues;
            var rightHandValues = rightHandPart.Values;

            var columnsCount = Math.Min(matrix.RowsCount, rightHandPart.Length);

            var elementIndex = begins[columnsCount] - 1;

            for (var columnIndex = columnsCount - 1; 0 < columnIndex; columnIndex--)
            {
                var rowStart = begins[columnIndex];

                var currentValue = rightHandValues[columnIndex] / diagonalValues[columnIndex];
                rightHandValues[columnIndex] = currentValue;

                while (rowStart <= elementIndex)
                {
                    var rowIndex = indexes[elementIndex];
                    rightHandValues[rowIndex] -= triangleValues[elementIndex] * currentValue;

                    elementIndex--;
                }
            }

            rightHandValues[0] /= diagonalValues[0];
        }

        /// <summary>
        ///     Solves a system of linear algebraic equations with the specified matrix and the right-hand part and
        ///     stores the result in the right-hand part vector.
        /// </summary>
        /// <param name="matrix">
        ///     The matrix of system.
        /// </param>
        /// <param name="rightHandPart">
        ///     The right-hand part vector as input and the result vector as output.
        /// </param>
        /// <exception cref="ArgumentNullException">
        ///     <paramref name="matrix"/> is <see lagnword="null"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="rightHandPart"/> is <see lagnword="null"/>.
        ///     </para>
        /// </exception>
        /// <remarks>
        ///     The size of system to solve is determined as minimal value from rows count of <paramref name="matrix"/> and
        ///     length of <paramref name="rightHandPart"/>.
        /// </remarks>
        public static void SolveWithLowTriangle(
            this CrcsSparseSymmetricMatrix<decimal> matrix, DenseVector<decimal> rightHandPart)
        {
#if SAFE

            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (rightHandPart == null)
            {
                throw new ArgumentNullException("rightHandPart");
            }

#endif

            var begins = matrix.Begins;
            var indexes = matrix.Indexes;
            var triangleValues = matrix.TriangleValues;
            var rightHandValues = rightHandPart.Values;

            var elementIndex = begins[1];

            var rowsCount = Math.Min(matrix.RowsCount, rightHandPart.Length);
            for (var rowIndex = 1; rowIndex < rowsCount; rowIndex++)
            {
                var rowEnd = begins[rowIndex + 1];

                var currentValue = rightHandValues[rowIndex];

                while (elementIndex < rowEnd)
                {
                    var columnIndex = indexes[elementIndex];
                    currentValue -= rightHandValues[columnIndex] * triangleValues[elementIndex];

                    elementIndex++;
                }

                rightHandValues[rowIndex] = currentValue;
            }
        }

        /// <summary>
        ///     Solves a system of linear algebraic equations with the specified matrix and the right-hand part and
        ///     stores the result in the right-hand part vector.
        /// </summary>
        /// <param name="matrix">
        ///     The matrix of system.
        /// </param>
        /// <param name="rightHandPart">
        ///     The right-hand part vector as input and the result vector as output.
        /// </param>
        /// <exception cref="ArgumentNullException">
        ///     <paramref name="matrix"/> is <see lagnword="null"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="rightHandPart"/> is <see lagnword="null"/>.
        ///     </para>
        /// </exception>
        /// <remarks>
        ///     The size of system to solve is determined as minimal value from rows count of <paramref name="matrix"/> and
        ///     length of <paramref name="rightHandPart"/>.
        /// </remarks>
        public static void SolveWithLowTriangleAndDiagonal(
            this CrcsSparseSymmetricMatrix<decimal> matrix, DenseVector<decimal> rightHandPart)
        {
#if SAFE

            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (rightHandPart == null)
            {
                throw new ArgumentNullException("rightHandPart");
            }

#endif

            var begins = matrix.Begins;
            var indexes = matrix.Indexes;
            var triangleValues = matrix.TriangleValues;
            var diagonalValues = matrix.DiagonalValues;
            var rightHandValues = rightHandPart.Values;

            var elementIndex = begins[1];

            rightHandValues[0] /= diagonalValues[0];

            var rowsCount = Math.Min(matrix.RowsCount, rightHandPart.Length);
            for (var rowIndex = 1; rowIndex < rowsCount; rowIndex++)
            {
                var rowEnd = begins[rowIndex + 1];

                var currentValue = rightHandValues[rowIndex];

                while (elementIndex < rowEnd)
                {
                    var columnIndex = indexes[elementIndex];
                    currentValue -= rightHandValues[columnIndex] * triangleValues[elementIndex];

                    elementIndex++;
                }

                rightHandValues[rowIndex] = currentValue / diagonalValues[rowIndex];
            }
        }

        /// <summary>
        ///     Subtracts the second matrix from the specified elements of the first matrix.
        /// </summary>
        /// <param name="minuend">
        ///     The minuend matrix.
        /// </param>
        /// <param name="minuendRowsIndexes">
        ///     The collection of zero-based indexes of rows of minuend elements.
        /// </param>
        /// <param name="minuendColumnsIndexes">
        ///     The collection of zero-based indexes of columns of minuend elements.
        /// </param>
        /// <param name="subtrahend">
        ///     The subtrahend matrix.
        /// </param>
        /// <exception cref="ArgumentNullException">
        ///     <paramref name="minuend"/> is <see langword="null"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="subtrahend"/> is <see langword="null"/>.
        ///     </para>
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="minuendRowsIndexes"/> is <see langword="null"/>.
        ///     </para>
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="minuendColumnsIndexes"/> is <see langword="null"/>.
        ///     </para>
        /// </exception>
        /// <exception cref="ElementIndexIsOutOfRangeException">
        ///     An element of <paramref name="minuendRowsIndexes"/> is negative or greater than or equal to
        ///     rows count of <paramref name="minuend"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         An element of <paramref name="minuendColumnsIndexes"/> is negative or greater than or equal to
        ///         columns count of <paramref name="minuend"/>.
        ///     </para>
        /// </exception>
        /// <exception cref="ReadOnlyMatrixElementException">
        ///     An element to assign specified by <paramref name="minuendRowsIndexes"/> and
        ///     <paramref name="minuendColumnsIndexes"/> is not in the portrait of <paramref name="minuend"/>.
        /// </exception>
        /// <remarks>
        ///     The number of rows of <paramref name="subtrahend"/> to subtract is determined as
        ///     min(<paramref name="minuendRowsIndexes"/>.Length, <paramref name="subtrahend"/>.RowsCount).
        ///     The number of columns of <paramref name="subtrahend"/> to subtract is determined as
        ///     min(<paramref name="minuendColumnsIndexes"/>.Length, <paramref name="subtrahend"/>.ColumnsCount).
        /// </remarks>
        public static void Subtract(
            this CrcsSparseSymmetricMatrix<decimal> minuend,
            int[] minuendRowsIndexes,
            int[] minuendColumnsIndexes,
            DenseSymmetricMatrix<decimal> subtrahend)
        {
#if SAFE

            if (minuend == null)
            {
                throw new ArgumentNullException("minuend");
            }

            if (minuendRowsIndexes == null)
            {
                throw new ArgumentNullException("minuendRowsIndexes");
            }

            if (minuendColumnsIndexes == null)
            {
                throw new ArgumentNullException("minuendColumnsIndexes");
            }

            if (subtrahend == null)
            {
                throw new ArgumentNullException("subtrahend");
            }

#endif

            var minuendSize = minuend.RowsCount;
            var begins = minuend.Begins;
            var indexes = minuend.Indexes;
            var minuendDiagonalValues = minuend.DiagonalValues;
            var minuendValues = minuend.TriangleValues;
            var subtrahendValues = subtrahend.Values;

            var subtrahendElementIndex = 0;

            var sizeToSubtract =
                Math.Min(subtrahend.RowsCount, Math.Min(minuendRowsIndexes.Length, minuendColumnsIndexes.Length));
            for (var rowIndex = 0; rowIndex < sizeToSubtract; rowIndex++)
            {
                var rowElementIndex = minuendRowsIndexes[rowIndex];

#if SAFE

                if ((rowElementIndex < 0) || (minuendSize <= rowElementIndex))
                {
                    var argumentName = string.Format(CultureInfo.CurrentCulture, "minuendRowsIndexes[{0}]", rowIndex);
                    throw new ElementIndexIsOutOfRangeException(argumentName, rowElementIndex, 0, minuendSize - 1);
                }

#endif

                var rowBegin = begins[rowElementIndex];
                var rowEnd = begins[rowElementIndex + 1];

                for (var columnIndex = 0; columnIndex <= rowIndex; columnIndex++)
                {
                    var columnElementIndex = minuendColumnsIndexes[columnIndex];

#if SAFE

                    if ((rowIndex == 0) && ((columnElementIndex < 0) || (minuendSize <= columnElementIndex)))
                    {
                        var argumentName = string.Format(CultureInfo.CurrentCulture, "minuendColumnsIndexes[{0}]", columnIndex);
                        throw new ElementIndexIsOutOfRangeException(argumentName, columnElementIndex, 0, minuendSize - 1);
                    }

                    if (rowElementIndex < columnElementIndex)
                    {
                        throw new ReadOnlyMatrixElementException("minuend", rowElementIndex, columnElementIndex);
                    }

#endif

                    if (rowElementIndex == columnElementIndex)
                    {
                        minuendDiagonalValues[rowElementIndex] -= subtrahendValues[subtrahendElementIndex];
                        subtrahendElementIndex++;
                        continue;
                    }

                    var minuendElementIndex = rowBegin;
                    while ((minuendElementIndex < rowEnd) && (indexes[minuendElementIndex] < columnElementIndex))
                    {
                        minuendElementIndex++;
                    }

                    if ((minuendElementIndex < rowEnd) && (indexes[minuendElementIndex] == columnElementIndex))
                    {
                        minuendValues[minuendElementIndex] -= subtrahendValues[subtrahendElementIndex];
                    }
                    else
                    {
#if SAFE

                        throw new ReadOnlyMatrixElementException("minuend", rowElementIndex, columnElementIndex);

#endif
                    }

                    subtrahendElementIndex++;
                }
            }
        }

        /// <summary>
        ///     Subtracts the second matrix multiplied by the factor from the specified elements of the first matrix.
        /// </summary>
        /// <param name="minuend">
        ///     The  minuend matrix.
        /// </param>
        /// <param name="minuendRowsIndexes">
        ///     The collection of zero-based indexes of rows of minuend elements.
        /// </param>
        /// <param name="minuendColumnsIndexes">
        ///     The collection of zero-based indexes of columns of minuend elements.
        /// </param>
        /// <param name="subtrahend">
        ///     The subtrahend matrix.
        /// </param>
        /// <param name="subtrahendFactor">
        ///     The factor of the subtrahend matrix.
        /// </param>
        /// <exception cref="ArgumentNullException">
        ///     <paramref name="minuend"/> is <see langword="null"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="subtrahend"/> is <see langword="null"/>.
        ///     </para>
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="minuendRowsIndexes"/> is <see langword="null"/>.
        ///     </para>
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         <paramref name="minuendColumnsIndexes"/> is <see langword="null"/>.
        ///     </para>
        /// </exception>
        /// <exception cref="ElementIndexIsOutOfRangeException">
        ///     An element of <paramref name="minuendRowsIndexes"/> is negative or greater than or equal to
        ///     rows count of <paramref name="minuend"/>.
        ///     <para>
        ///         --or--
        ///     </para>
        ///     <para>
        ///         An element of <paramref name="minuendColumnsIndexes"/> is negative or greater than or equal to
        ///         columns count of <paramref name="minuend"/>.
        ///     </para>
        /// </exception>
        /// <exception cref="ReadOnlyMatrixElementException">
        ///     An element to assign specified by <paramref name="minuendRowsIndexes"/> and
        ///     <paramref name="minuendColumnsIndexes"/> is not in the portrait of <paramref name="minuend"/>.
        /// </exception>
        /// <remarks>
        ///     The number of rows of <paramref name="subtrahend"/> to subtract is determined as
        ///     min(<paramref name="minuendRowsIndexes"/>.Length, <paramref name="subtrahend"/>.RowsCount).
        ///     The number of columns of <paramref name="subtrahend"/> to subtract is determined as
        ///     min(<paramref name="minuendColumnsIndexes"/>.Length, <paramref name="subtrahend"/>.ColumnsCount).
        /// </remarks>
        public static void Subtract(
            this CrcsSparseSymmetricMatrix<decimal> minuend,
            int[] minuendRowsIndexes,
            int[] minuendColumnsIndexes,
            DenseSymmetricMatrix<decimal> subtrahend,
            decimal subtrahendFactor)
        {
#if SAFE

            if (minuend == null)
            {
                throw new ArgumentNullException("minuend");
            }

            if (minuendRowsIndexes == null)
            {
                throw new ArgumentNullException("minuendRowsIndexes");
            }

            if (minuendColumnsIndexes == null)
            {
                throw new ArgumentNullException("minuendColumnsIndexes");
            }

            if (subtrahend == null)
            {
                throw new ArgumentNullException("subtrahend");
            }

#endif

            var minuendSize = minuend.RowsCount;
            var begins = minuend.Begins;
            var indexes = minuend.Indexes;
            var minuendDiagonalValues = minuend.DiagonalValues;
            var minuendValues = minuend.TriangleValues;
            var subtrahendValues = subtrahend.Values;

            var subtrahendElementIndex = 0;

            var sizeToSubtract =
                Math.Min(subtrahend.RowsCount, Math.Min(minuendRowsIndexes.Length, minuendColumnsIndexes.Length));
            for (var rowIndex = 0; rowIndex < sizeToSubtract; rowIndex++)
            {
                var rowElementIndex = minuendRowsIndexes[rowIndex];

#if SAFE

                if ((rowElementIndex < 0) || (minuendSize <= rowElementIndex))
                {
                    var argumentName = string.Format(CultureInfo.CurrentCulture, "minuendRowsIndexes[{0}]", rowIndex);
                    throw new ElementIndexIsOutOfRangeException(argumentName, rowElementIndex, 0, minuendSize - 1);
                }

#endif

                var rowBegin = begins[rowElementIndex];
                var rowEnd = begins[rowElementIndex + 1];

                for (var columnIndex = 0; columnIndex <= rowIndex; columnIndex++)
                {
                    var columnElementIndex = minuendColumnsIndexes[columnIndex];

#if SAFE

                    if ((rowIndex == 0) && ((columnElementIndex < 0) || (minuendSize <= columnElementIndex)))
                    {
                        var argumentName = string.Format(CultureInfo.CurrentCulture, "minuendColumnsIndexes[{0}]", columnIndex);
                        throw new ElementIndexIsOutOfRangeException(argumentName, columnElementIndex, 0, minuendSize - 1);
                    }

                    if (rowElementIndex < columnElementIndex)
                    {
                        throw new ReadOnlyMatrixElementException("minuend", rowElementIndex, columnElementIndex);
                    }

#endif

                    if (rowElementIndex == columnElementIndex)
                    {
                        minuendDiagonalValues[rowElementIndex] -= subtrahendFactor * subtrahendValues[subtrahendElementIndex];
                        subtrahendElementIndex++;
                        continue;
                    }

                    var minuendElementIndex = rowBegin;
                    while ((minuendElementIndex < rowEnd) && (indexes[minuendElementIndex] < columnElementIndex))
                    {
                        minuendElementIndex++;
                    }

                    if ((minuendElementIndex < rowEnd) && (indexes[minuendElementIndex] == columnElementIndex))
                    {
                        minuendValues[minuendElementIndex] -= subtrahendFactor * subtrahendValues[subtrahendElementIndex];
                    }
                    else
                    {
#if SAFE

                        throw new ReadOnlyMatrixElementException("minuend", rowElementIndex, columnElementIndex);

#endif
                    }

                    subtrahendElementIndex++;
                }
            }
        }

        #endregion

        #endregion
    }
}